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Abstract 

In the reconstruction process of unknown multiple scattering objects in inverse 
medium scattering problems, the first important step is to effectively locate some 
approximate domains that contain all inhomogeneous media. Without such an 
effective step, one may have to take a much larger computational domain than ac- 
tually needed in the reconstruction of all scattering objects, thus resulting in a huge 
additional computational efforts. In this work we propose a simple and efficient 
multilevel reconstruction algorithm to help locate an accurate position and shape 
of each inhomogeneous medium. Then other existing effective but computation- 
ally more demanding reconstruction algorithms may be applied in these initially 
located computational domains to achieve more accurate shapes of the scatter and 
the contrast values over each medium domain. The new algorithm exhibits several 
strengths: robustness against noise, requiring less incidences, fast convergence, flex- 
ibility to deal with scatterers of special shapes, and advantages in computational 
complexity. 
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1 Introduction 

This works investigates the numerical identification of inhomogeneous medium scat- 
terers by scattered fields. The inverse scattering problem can find wide applications 
in medicine, geophysics, biological studies. A large variety of numerical reconstruction 
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methods are available in literature, such as the time-reversal multiple signal classifica- 
tion (MUSIC) method [HIEI, the contrast source inversion (CSI) method [HHTTJII], the 
continuation method [2], the subspace-based optimization method (SBOM) |11|5], the hn- 
ear sampling or probing methods (LSM) [HI HSl [H], the parallel radial bisection method 
(PRBM) [13], etc. In order to carry out any of these methods for the reconstruction of 
unknown multiple scattering objects, the first important step is to effectively locate some 
approximate domains that contain all scattering objects. Without such an effective step, 
one may have to take a much larger computational domain than actually needed in the 
reconstruction of all scattering objects. In particular, when multiple separated objects are 
present, and at least two of them are far away from each other, then one may need to set 
an initial sampling or computational domain to be much larger in order to ensure a safe 
covering of all such objects. One may easily have an initial computational domain with 
an area or a volume of 30 or 40 times as large as the actual region required to cover all 
unknown inhomogeneous medium objects. A much larger computational domain results 
usually in a huge additional computational effort for the entire numerical reconstruction 
process, considering the severe ill-posedness and strong nonlinearity of inverse medium 
scattering problems. 

So it is of great significance for the reconstruction process of an inverse medium prob- 
lem to have an effective step that helps locate the initial regions that cover each of the 
scattering object. In addition, this first step should be less expensive computationally 
and easy to implement numerically. It is mostly challenging to realize this target, and to 
provide an acceptable initial location of each scattering object at the same time. A direct 
sampling method was proposed recently in [7] for the purpose. The algorithm is computa- 
tionally very cheap as it involves computing only the inner product of the scattered field 
with fundamental solutions located at sampling points. In this work, we will propose a 
new algorithm for the purpose, and it is completely different from the one in [7] . This new 
algorithm is an iterative one, also very cheap; only three matrix-vector multiplications 
are needed at each iteration, without any matrix inversion or solutions of linear systems 
involved. Most interestingly, the algorithm can first separate all disjoint objects quickly, 
usually in a few iterations, then refine its approximation successively and finally provide 
a good approximate domain for each separate object. 

It is worth mentioning that the multilevel algorithm to be presented here is essentially 
different in nature from the multilevel linear sampling method developed in [11]: the 
new method is much less sensitive to the so-called cut-off values, it works with much 
less incident fields, and it does not need to solve ill-posed far-field equations at sampling 
points. In addition, the new algorithm is robust in the presence of noise and less sensitive 
to the noise level. More importantly, unlike most existing methods, the new method does 
not involve any optimization process or matrix inversions, so it can be viewed as a direct 
sampling method. Another nice feature of the new algorithm is that it is self-adaptive, 
that is, at each iteration it can remedy the possible errors from the previous iterations. 
With an effective initial location of each object, we may then apply any existing efficient 
but computationally more demanding methods, e.g., the methods in [2] [S] [T7], for further 
refinement of the estimated shape of each scattering object as well as for recovery of the 
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contrast profiles of different media. 



2 Problem description 

Consider an inverse scattering problem where the scatterer Q, possibly consisting of 
several separated disjoint components, is located in a homogeneous background medium 

{d = 2,3). We assume that the scattered obstacles are illuminated successively by 
a number of plane wave incident fields u^'^lx), j = 1,2, ■■ ■ ,iVj. For each plane wave 
incidence, the scattered field u''"^°'{Xg) is measured by the receivers at locations xl, ■ ■ ■ , 
x%- ; see Figure [H for the incidences and receivers located on a circle S. 




Figure 1 : Geometrical model of the scattering problem. 

The inverse scattering problem is to determine the contrast function or index of refraction, 
x{x) for any point x varying in the scatterer Q, given a set of NiNg scattering data m'^™. 
The contrast x has a very important property, namely it vanishes outside the scattering 
obstacles. For each incident field, the total field Uj satisfies the Helmholtz equation [6]: 

Auj{x) + k'^ixix) + l)uj{x) = 0, X e W^, (2.1) 

where k is the wavenumber of the homogeneous background media. The total field Uj can 
be represented by the intergral equation [6], 

u^{x)=uf^{x) + k'^ I g{x,x')x{x')uj{x')dv{x'), (2.2) 

where g{x,x') is the Green's function of the homogeneous background medium: 

, lH;,'\k\x-x'\) for d = 2, 

where ifp^^ is the zero-order Hankel function of first kind. We note that the total field u 
may stand for the acoustic pressure in an acoustic scattering problem, or for the electric 
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field vector in an electromagnetic scattering, or for tlie particle-velocity vector in an 
elastodynamic scattering. Tlie scattered field is measured on the boundary of a domain 
S, which is sitting outside the scatterer Q. As the scatterer Q is unknown, we shall 
introduce a sampling domain D that completely cover the scatterer Q. As the contrast 
function x vanishes outside Q, with the help of f l2.2p we can write the scattered field as 

up{x) = Uj{x) - uf\x) = k^ [ gix, x')xix')uj{x')dv{x'), x e S. (2.3) 

For sake of convenience, we shall often introduce the contrast source function 

Wjix) = xix)ujix), xeD. (2.4) 
Then we can write (12. 2 p and (12. 3p in the following more compact forms 

wj{x) = x{x)uf'{x) + x{x){GDWj){x), X G D, (2.5) 

and 



uT^ix) = {GsWj){x), xeS, (2.6) 
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where Gd and Gs are two integral operators given by 

{Gdw){x) = k'^ I g{x,x')w{x')dv{x') Va;GD, 
Jd 

iGsw){x) = k'^ [ g{x,x')w{x')dv{x') Vx G S. 
Jd 

Equations (12. 5p and (12. 6p will be two fundamental equations for our proposed multilevel 
initialization algorithm. 

3 Approximate contrast source by backpropagation 

We can easily see that the support of the contrast source function w = x^ describes 
the exact locations and geometries of all the inhomogeneous media, which generate a 
scattered field m*'^". The aim of this work is to propose a fast and less expensive algo- 
rithm that can help locate all the inhomogeneous media and provide good initial guesses 
for some computationally more demanding iterative algorithms to find more accurate 
approximations of the contrast function x- 

Our algorithm will rely on the approximate contrast source obtained by backpropa- 
gation. Backprogation is widely used in inverse medium scatterings, see [T7] [TU] and the 
references therein. In this section we shall give a rigorous mathematical explanation of 
the approximate contrast source by backpropagation. Let {■,-)l'^{s) and {■,-)l'^{d) be the 
scalar products respectively in L^{S) and L'^{D), and G*:L'^{S) L^{D) be the adjoint 
of operator G^: L^{D) — > L^{S). G* is called the backpropagation operator and given by 

iG*sw){x) = k'^ [ g{x,x')w{x') rfs(x') Vx G D. 
Js 
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We shall need the following backpropagation subspace of L^{D), 

Vb = span {G*y^''} , 

which is formed by all the fields generated by the backpropagation G* on the scattered 
data u^^"'. It follows from fl2.6p that 

u'^^ix) = {Gsw){x), xeS. (3.1) 

The backpropagation is to seek a best approximate solution Wb to the equation ( 13. ip in 
the backpropagation subspace Vf,, namely 

\\u''''-GsWb\\l2^s) = min \\u'''''' - GsVbWl^sy (3.2) 

It is easy to see that the solution Wb to ( 13. 2p solves the variational system: 

{u'''''-GsWb,GsVb)L^^s) = yvbeVb, (3.3) 

or equivalently, 

{GsWb, GsVb)L^s) = {Glu'"", Vb)L^(D) ^VbEVb. (3.4) 
As Wb,Vb G Vb, we can write 

Wb = X Gy^" , Vb = fi G*y' (3.5) 
for some constants X, fi. Substituting the two expressions into (13. 4p we obtain 

1 1 sea 1 1 2 

which gives the approximate contrast source by backpropagation: 

I I I I 2 

-'= 110.0: J G-y-. (3.7) 



4 A multilevel algorithm 

In this section we propose a fast multilevel algorithm to find the locations and geome- 
tries of all the inhomogeneous media, which are described by the contrast function x 
(12. ip . The algorithm proceeds iteratively, and carries out two important steps at each 
iteration based on the two fundamental equations (12. 5 p and (12. 6p . namely the state and 
field equations. In the first step, we apply the backpropagation technique to compute an 
approximate contrast source Wj corresponding to each incident u^^'^ (j = 1, 2, ■ ■ ■ , TVj). It 
follows from (13. 7p that this approximation is given by 

1 1 ySca 1 1 2 

^.= 11^ G*sur, j = l,2,...,iV.. (4.1) 
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With these approximate contributions Wj of the exact contrast source w corresponding to 
each incident m*"'^, we approximate the contrast x pointwise by minimizing the residual 
equation corresponding to the state equation fl2.5l) . namely 



which yields an explicit formula to compute an approximate contrast value x(^) every 
point X E D when an approximate contrast source Wj is available: 



where the overbear denotes the complex conjugate. 

Clearly both (14. ip and (14.31) are rather crude mostly, and may provide rather poor 
approximations for the exact contrast source w and contrast profile x- But as it will 
be seen, when we combine these two poor approximations in a novel manner with some 
multilevel technique, it generates a very efficient and robust algorithm for locating an 
accurate position and shape of each inhomogeneous medium. 

Our basic idea is also based on another simple observation. We know that the exact 
contrast function vanishes outside the scatterer f2, so its support provides the lo- 
cation and shape of the scatterer f2, which is formed by all the inhomogeneous media. 
This observation, along with the previous two explicit evaluation formulae (14. ip and (14. 3 p 
and a novel multilevel technique, forms the foundation of our new multilevel sampling 
algorithm. 

For the description of the algorithm, we first introduce two new concepts, the smallest 
distance and the first gap interval with index M. For a given finite positive sequence, 
{ai, a2, ■ ■ ■ , Cm}, its smallest distance is the positive smallest one among all the distances 
between two neighboring elements, namely dist{ai,ai^i), z = 1, 2, ■ ■ ■ , m — 1. Among all 
these m — 1 distances, if there exists some j such that 2 < j < m — 1 and the distance 
dist{aj, ctj+i) is M times larger than the smallest distance of the sequence {ai, 02, ■ ■ ■ , aj}, 
then [aj, aj+i] is called a gap interval. The first such interval is called the first gap interval. 

Now we are ready to state our new algorithm. 
Multilevel Sampling Algorithm. 

1. Choose a sampling domain D that contains the scatterer Q; 

Select a uniform (coarse) mesh on D, consisting of rectangular (2D) or cubic (3D) 
elements; write the mesh as Di, 

Select a tolerance e and an index M; set an initial cut-off value cq := and k := 1. 

2. Compute an approximate value of the contrast Xfc(^) ^it each grid point x G Dk, 
using the formulae (14. ip and (14. 3p . Then do the following: 




(4.2) 




(4.3) 
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2.1 Order all the values of Xk{x) satisfying Xkix) > Ck-i into a non- decreasing 
sequence; 

Find the first gap interval of the sequence with index M; 

Choose the right endpoint of this first gap interval as the next cut-off value c^. 

2.2 If Xk{x) > Cfc at a grid point 03, select all the grid points of the elements which 
share x as one of its vertices; 

Remove all the grid points in Z)^, which are not selected; 
Update Dk by all those selected grid points. 

3. If |cfe — Cfc-i| < e, go to Step 4; 

otherwise refine the mesh Dk to get -Dfc+i; set k := k + 1 and go to Step 2. 

4. Output all grid points in for the supports of all inhomogeneous media; 
Output the approximate contrast values Xk{'^) at all grid points of Dj.. 

We can easily see that the Multilevel Sampling Algorithm does not involve any opti- 
mization process or matrix inversions, and its major cost is to update the contrast values 
using the explicit formulae (14. ip and (14. 3 p at each iteration, and the computational sam- 
pling domain Dk shrinks as the iteration goes. So the algorithm is rather simple and 
less expensive. In addition, as the cut-off values are basically to distinguish the homo- 
geneous background medium where x(^) vanishes and the inhomogeneous media where 
should be essentially different from 0, so our cut-off values are rather easy to choose 
and insensitive to the size and physical features of scatterers. In fact, the cut-off value 
can start simply with zero, then it is updated automatically with the iteration. As we 
shall see from numerical examples in the next section, the algorithm works well with few 
incidents, even with one; and it is self-adaptive, namely it can recover some elements that 
have been removed at the previous iterations due to the computational errors. In terms 
of these aspects the new Multilevel Sampling Algorithm outperforms the popular linear 
sampling methods [15], including the improved multilevel variant [TT] . 

Discretization. We end this section with a brief discussion about numerical dis- 
cretization of the integrations involved in the above multilevel algorithm. We illustrate 
only the discretization of (12. 5p and (12. 6p . as all other integrations involved in the algorithm 
can be approximated similarly. To do so, we divide the domain D into smaller rectangular 
or cubic elements, whose centers are denoted as a^i, 0^2, . . . ,xm- Using the coupled-dipole 
method (CDM) or discrete dipole approximation (DDA) [3l[9], we can discretize (12. 5p by 

Wj{Xm) = xiXm)uf{Xm) + k'^Xi^m) AnQ^Xm, Xn)Wj{Xn), m = 1, 2, . . . , M, (4.4) 

where An is the area or volume of the n-th element. Similarly, we can discretize equation 
(12. 6 P at every point a; G S by 

AI 

^^^(a;) = A;2^A„(7(x,a;„K(ccO for j = 1, 2, ■ ■ ■ , A^^. (4.5) 

n=l 



7 



5 Numerical simulations 



In this section we present several examples to test the effectiveness of our multilevel 
sampling algorithm. We first list the parameters that are used in our numerical simu- 
lations: the wave number k and wave length A are taken to be = 27r and A = 1; the 
number of incidences and receivers are set to be A^j = 6 and Ng = 30 respectively, both 
equally distributed on the circle of radius 5A; the index M of the first gap interval and 
the tolerance parameter e are chosen to be 100 and 10"'^ respectively. In all the numerical 
simulations, random noises are added to the exact scattering data in the following form: 

u^^'^ix) := uf''{x)[l + + ^r2A^))], J = 1, 2, ■ ■ ■ , iV^ 

where rij{x) and r2j{x) are two random numbers varying between -1 and 1, and ^ 
corresponds to the level of the noise, which is usually taken to be 10% unless specified 
otherwise. All the programs in our experiments are written in MATLAB and run on a 
2.83GHz PC with 4GB memory. 

5.1 Two-dimensional reconstructions 

Example 1. This example shows a scatterer Q consisting of two squares of side length 
0.3A, located respectively at (— 0.3A, — 0.3A) and (0.3A,0.3A), with their contrast values 
being 1 and 2 respectively; see the two red squares in Figure [2]^a). We take the sampling 
domain D = [— 1.2A, 1.2A] x [— 1.2A, 1.2A], which is quite large compared to the scatterer 
fl, with its area of 64 times of the area of one scatterer component. More importantly, we 
see that these two small objects are quite close to each other. The mesh refinement during 
the multilevel algorithm is carried out based on the rule hk = 0.4A/2'^, where k is the 
k-th refinement, and and are respectively the mesh sizes of the initial mesh and the 
mesh after the kth refinement. The numerical reconstructions are shown in Figure [2|^b) - 
Wid) respectively for the 1st, 3rd and 5th iteration. One can observe from the figures 
that the algorithm converges very fast and provides accurate locations of the two medium 
components in only 5 iterations. Moreover, we can see an important advantage of the 
algorithm, i.e., it can separate the disjoint medium components quickly. 

Example 2. This example is the same as Example 1, except that the contrast values 
of two medium components are now variable functions, namely 

, , . 7r(10|x| - 1.5) . 7i(10\y\ - 1.5) 
X[x,y) = sm ^ sm ^ . 

The numerical reconstructions are shown in Figure [3t^b)- [3](d) for the 1st, 3rd and 6th 
iteration. Again, we observe from the figures that the algorithm converges fast, provides 
very satisfactory locations of the two medium components in only 6 iterations, and it can 
separate the disjoint medium components quickly. 

Example 3. This example considers a scatterer of a thin annulus with the inner 
and outer radii being 0.3A and 0.5A respectively and centered at the origin. The contrast 
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(a) 



(b) 



(c) 



(d) 



Figure 2: (a) The initial (coarse) mesh on the sampling domain; (h)-(d) Reconstructions 
at the 1st, 3rd and 5th iteration. 



value x{'^) is 2 inside the thin annulus. The samphng domain D is taken to be a square 
of side length with 5.6A, as shown in Figure Hl^a). 

It is easy to see the sampling domain D has an area of about 62 times as large as 
the one of the annulus, and the annulus has a very thin thickness, i.e., 0.2A. The mesh 
refinement during the multilevel algorithm is carried out based on the rule hk = 0.4A/2'^, 
where k indicates the k-th refinement and hk is the mesh size after the kth refinement. The 
numerical reconstructions are shown in Figure |4]^b)-|ll^d) for the 1st, 3rd and 4th iteration. 
Same as for the previous two examples, the reconstructions are quite satisfactory and the 
accurate locations of the scatterer can be achieved. 

5.2 Reconstruction for the contrast function x 

There are many numerical reconstruction methods available in the literature for the 
contrast profile function x, that are robust and more accurate than our new multilevel 
method. But these methods are usually more complicated technically and much more 
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(a) 



(b) 



(c) 



(d) 



Figure 3: (a) The first coarse mesh on the sampling region; (h)-(d) Reconstructions at the 
1st, 3rd and 6th iteration. 



demanding computationally, as they mostly involve nonlinear optimizations and matrix 
inversions. Without a reasonably good initial location for each inhomogeneous medium, 
we may have to take a much larger sampling domain than required so these methods can 
be extremely time consuming, especially in three dimensions. Using the newly proposed 
multilevel algorithm in Section HI we can first locate a much smaller sampling domain 
than usual (or the one we originally selected) in a numerical reconstruction for the con- 
trast X- Then we can apply any existing reconstruction algorithms for more accurate 
reconstructions, starting with an initial sampling domain provided by the new multilevel 
algorithm. This may save us a great fraction of the entire computational costs. 

Next we show numerical tests by combining our new multilevel method with the 
popular extended contrast source inversion (ECSI) method [17]. Firstly, we consider 
the same scatterer VL and the set-ups as in Example 1 of Section I5.lt see Figure [5](a). 
Then we apply the ECSI method [17] with mesh size h = O.OIA to the reconstructed 
domain (cf. Figure [2]^d)) by the multilevel algorithm. The reconstruction result is shown 
in Figure [5]^b). Clearly, one can see that both the location and the values of the contrast 
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(a) 



(b) 




(c) 



(d) 



Figure 4: (a) The first coarse mesh on the sampling region; (h)-(d) Reconstructions at 
the 1st, 3rd and 4th iteration. 



are well reconstructed. 

Secondly, we consider the same scatterer Q and the set-ups as in Example 3 of Section 
I5.1j see Figure (E])(c). We apply the ECSI method [T7] with mesh size h = 0.025A to 
the reconstructed domain (cf. Figure 111(d)) obtained by the multilevel algorithm. The 
reconstruction result is shown in Figure El^d). Clearly, both of the location and the values 
of the contrast are well reconstructed. 



5.3 Three-dimensional reconstructions 

Example 4. This example tests a three-dimensional scatterer Q consisting of two 
small cubic components: 

ni = [-0A5\, -0.15A]^ 1^2 = [0.15A, 0.45A]^ 

The two squares are quite close to each other, both with constant contrast values 2; see 
Figure ini^a). We take the sampling domain to be D = [— 1.2A, 1.2A]^, which is about 500 
times of the volume of Qi or ^22- 
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(c) (d) 

Figure 5: (a) True scatterer in Example 1; (b) Reconstruction by combining the multilevel 
method with ECSI for Example 1; (c) True scatterer in Example 3; (d) Reconstruction by 
combining the multilevel method with ECSI for Example 3 



We take an initial mesh size of Hq = 0.8A in the multilevel algorithm. The mesh refine- 
ment during the multilevel algorithm is carried out based on the rule hk = 0.8A/2'^, where 
k is the k-th refinement, and hk is the mesh size after the kth refinement. The numerical 
reconstructions are shown in Figure [71 Same as for the previous two-dimensional exam- 
ples, the reconstructions are quite satisfactory and the accurate locations of the scatterer 
can be achieved, and two inhomogeneous medium objects can be quickly separated. 

Example 5. In this test we consider a torus scatterer (see Figure Et^b)), with a 
contrast value 2. The torus has the following representation, 

where r = O.IA and R = 0.4A {R is the radius from the center of the hole to the center 
of the torus tube, r is the radius of the tube). The sampling doman is taken to be 
D = [— 1.2A, 1.2A]^, which is about 170 times of the volume of the torus. 
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(b) 



(a) 



Figure 6: Scatterers imbedded in a large sampling domain: two cubic components close 
to each other (a); a torus (b) 




(c) (d) 
Figure 7: Numerical reconstructions by the first 4 iterations 
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(a) 



(b) 




(c) 



(d) 



Figure 8: Numerical reconstructions by the first 3 iterations; (d) is the same as (c), hut 
viewed in a different angle. 

We take an initial mesh size of Hq = 0.4A in the multilevel algorithm. The mesh 
refinement during the multilevel algorithm is carried out based on the rule hk = 0.4A/2'^, 
where k is the k-th refinement, and hk is the mesh size after the kth refinement. The 
numerical reconstructions are shown in Figure El Same as for the previous two-dimensional 
examples, the reconstructions are quite satisfactory and the accurate locations of the 
scatterer can be achieved. 

6 Concluding remarks 

This work proposes a multilevel sampling algorithm which helps locate an initial com- 
putational domain for the numerical reconstruction of inhomogeneous media in inverse 
medium scatterings. The algorithm is an iterative process which starts with a large sam- 
pling domain, and reduces the size of the domain iteratively based on the cut-off values, 
which are computed adaptively by using the updated contrast source strengths and con- 
trast values at each iteration. The iterative algorithm can be viewed actually as a direct 
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method, since it involves only matrix-vector operations and does not need any optimiza- 
tion process or to solve any large-scale ill-posed linear systems. The algorithm works with 
very few incident fields and its cut-off values are easy to compute and insensitive to the 
sizes and shapes of the scatterers, as well as the noise in the data. This is a clear ad- 
vantage of the algorithm over some popular existing sampling methods such as the linear 
sampling type methods, where the cut-off values are sensitive to the noise and difficult 
to choose, and the number of incident fields can not be small. In addition, the multilevel 
algorithm converges fast and can easily separate multiple disjoint scattering components, 
often with just a few iterations to find a satisfactory initial location of each object. An- 
other nice feature of the new algorithm is that it is self-adaptive, that is, it can remedy 
the possible errors from the previous levels at each current level. With an effective initial 
location of each object, we may then apply any existing efficient but computationally 
more demanding methods for further refinement of the estimated shape of each scattering 
object as well as for recovery of the contrast profiles of different media. 
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